#!/usr/bin/env python
#This tool allow users to plot SVM-prob ROC curve from data
from svmutil import *
from sys import argv, platform
from os import path, popen
from random import randrange , seed
from operator import itemgetter
from time import sleep
from math import exp
import copy

def myread(model_vector):
        prob_y = []
        prob_x = []
        for vector in model_vector:
             xi = {}
             value = vector[1]
             features = vector[0]
             count = 1 ; 
             for feature in features:
                 xi[count] = feature
                 count += 1 
             prob_y += [float(value)]
             prob_x += [xi]
        return (prob_y,prob_x)  


def get_pos_deci(train_y, train_x, test_y, test_x, param):
	model = svm_train(train_y, train_x, param)
	labels = model.get_labels()
	py, evals, deci = svm_predict(test_y, test_x, model)
	deci = [labels[0]*val[0] for val in deci]
	return deci,model,py


def get_cv_deci(prob_y, prob_x, param, nr_fold):
	if nr_fold == 1 or nr_fold==0:
		deci,model,py = get_pos_deci(prob_y, prob_x, prob_y, prob_x, param)
		return deci
	deci, model , py = [], [], []
	prob_l = len(prob_y)

	for i in range(prob_l):
		j = randrange(i,prob_l)
		prob_x[i], prob_x[j] = prob_x[j], prob_x[i]
		prob_y[i], prob_y[j] = prob_y[j], prob_y[i]

	for i in range(nr_fold):
		begin = i * prob_l // nr_fold
		end = (i + 1) * prob_l // nr_fold
		train_x = prob_x[:begin] + prob_x[end:]
		train_y = prob_y[:begin] + prob_y[end:]
		test_x = prob_x[begin:end]
		test_y = prob_y[begin:end]
		subdeci,submdel,subpy = get_pos_deci(train_y, train_x, test_y, test_x, param)
		deci += subdeci
                py += subpy
	return deci,py



def plot_roc(deci, label):
	#count of postive and negative labels
	db = []
	pos, neg = 0, 0 		
	for i in range(len(label)):
		if label[i]>0:
			pos+=1
		else:	
			neg+=1
		db.append([deci[i], label[i]])

	#sorting by decision value
	db = sorted(db, key=itemgetter(0), reverse=True)

	#calculate ROC 
	xy_arr = []
	tp, fp = 0., 0.			#assure float division
	for i in range(len(db)):
		if db[i][1]>0:		#positive
			tp+=1
		else:
			fp+=1
		xy_arr.append([fp/neg,tp/pos])

	#area under curve
	aoc = 0.			
	prev_x = 0
	for x,y in xy_arr:
		if x != prev_x:
			aoc += (x - prev_x) * y
			prev_x = x

        return aoc  	


def split_list(alist,w_p=1):
   length = len(alist)
   return [alist[i*length // w_p: (i+1)*length // w_p ] for i in range(w_p)]


def match(a,b):
    m = [ i for i , j in zip(a,b) if i==j ]
    r = len(m)/float((len(a)-len(m)))
    return r

def main():
         
         # Input file which contains info about the weight
         # weight_file = open('weight','r')
         # count = 0
         # val_w = []
         # pos_samples = 0
         # neg_samples = 0
         # for line in weight_file:
         #    val_w.append(int(line))
         #    if(val_w[count] == 1):
         #         pos_samples = pos_samples + 1
         #    elif(val_w[count] == -1):
         #         neg_samples = neg_samples + 1
         #    else:
         #         print "Value other than +- 1 Not Allowed"
         #         raise SystemExit
         #    count = count + 1
         # weight_file.close()           

         # total_samples = pos_samples + neg_samples
         # print "positive_samples = ",pos_samples
         # print "negative_samples = ",neg_samples 
         # weight_pos = 1./pos_samples
         # weight_neg = 1./neg_samples
         
         # print weight_pos , weight_neg
         # pos_weight = []
         # neg_weight = [] 
         # for i in range(total_samples):
         #    if(val_w[i] == 1):
         #         pos_weight.append(weight_pos)
         #    elif(val_w[i] == -1):
         #         neg_weight.append(weight_neg)  

         weight_file = open('weight','r')

         #### Number of Positive Samples ####
         pos_samples = 8903

         #### Number of Negative Samples ####
         neg_samples = 11706

         val_w = []

         total_samples = pos_samples + neg_samples

         for line in weight_file:
            val_w.append(float(line))

         pos_weight = val_w[0:pos_samples]

         neg_weight = val_w[pos_samples:total_samples]
         
         ow = []

         ow[0:pos_samples] = [1]*pos_samples
         ow[pos_samples:total_samples] = [-1]*neg_samples       
 
         # Number of Iterations
         T = 1

         # Number of features DIMENSION 
         d = 32
         
         # Number of local patches 
         lp = 435

         # Number of samples of each + and - to be taken for boosting  
         num_of_samples = 20*d

         print ">>>BOOSTING STARTED<<<"
         # BOOSTING WILL GO ON FROM HERE #

         optimum = []   
         f = open('optimum','r')
         for line in f:
             sval = int(line) 
             optimum.append(sval)
         f.close()

         for q in range(T):
           
           print "STAGE "+str(q)

           pdatabase = []
           ndatabase = []

         # Got Index after Sorting according to Weights    
           dpindex = [ i[0] for i in sorted(enumerate(pos_weight),key=lambda x:x[1]) ] 
           dnindex = [ i[0] for i in sorted(enumerate(neg_weight),key=lambda x:x[1]) ]

         # Initially take %num_of_samples%  Positive as well Negative Images   
           for i in dpindex[:num_of_samples]:
             fname = "pdata/"+str(i)
             input_file = open(fname,'r')
             for line in input_file:
                featureList =  line.split() 
                value = featureList[0]
                value = int(value)
                features = featureList[1:]
                fvalue = []   
                for feature in features:
                   fvalue.append(float(feature))
                pdatabase.append([fvalue,value]) 
             input_file.close()  
     
           for i in dnindex[:num_of_samples]:
             fname = "ndata/"+str(i)
             input_file = open(fname,'r')
             for line in input_file:
                featureList =  line.split() 
                value = featureList[0]
                value = int(value)
                features = featureList[1:]
                fvalue = []   
                for feature in features:
                   fvalue.append(float(feature))
                ndatabase.append([fvalue,value])  
             input_file.close() 

           pdbase = split_list(pdatabase,num_of_samples)
           ndbase = split_list(ndatabase,num_of_samples) 

           p_sep_data = []
           n_sep_data = []
           
 
           # Separate each Window Dimension   
           for i in range(lp):
             p_sep_data.append( [ a[i] for a in pdbase ] )
         
           for i in range(lp):
             n_sep_data.append( [ a[i] for a in ndbase ] )

 	   ############################################################
           ######## deepcopy makes entire new copy of the data ########
           ############################################################
      
 	   f_n_sep_data=copy.deepcopy(n_sep_data)
	   f_p_sep_data=copy.deepcopy(p_sep_data)

           #f_n_sep_data = n_sep_data
           #f_p_sep_data = p_sep_data
   

           for data in f_p_sep_data:
               t_count = 0
               for subwin in data:
                      for i in optimum:
                         subwin[0] = subwin[0] + p_sep_data[i][t_count][0]
                      t_count = t_count + 1
       
       
           for data in f_n_sep_data:
               t_count = 0 
               for subwin in data:
                      for i in optimum:
                         subwin[0] = subwin[0] + n_sep_data[i][t_count][0]
                      t_count = t_count + 1
 

           sep_data = [] 
           for i in range(lp):
             sep_data.append(f_p_sep_data[i] + f_n_sep_data[i])
          


           seed(0)
           
           criteria = 0 

           count = 0

           if criteria == 0: 
               max_auc = 0.0               
           else:
               max_index = 0 

           for vector in sep_data:
             if count not in optimum:
               prob_y,prob_x = myread(vector)
               fold = 5
	       param = ""
               org = [val[1] for val in vector]
               deci,py = get_cv_deci(prob_y,prob_x,param,fold)
               auc = plot_roc(deci, prob_y)
               
               if criteria == 0:
               	print "count = ",count," AUC = " , auc
               	if auc > max_auc:
                    max_auc = auc
                    max_index = count
                    max_py = py
               else:
                me = match(py,ow)
                print 'count = '+str(count)+'RATE = '+str(me)
                if me > max_rate:
                    max_rate = me 
                    max_index = count
                    max_py = py  
             count = count + 1
         
           if criteria == 0:    
           	print "MAX_AUC is ",max_auc
           else:
           	print "MAX_RATE is ",max_rate
     
           print "MAX_INDEX ",max_index

         ###################################################################################### 
         ##### Reweight the samples with respect to best classifier  found in the stage   #####
         ######################################################################################


         ###### positive samples ###### 
           for i in range(num_of_samples):
                 v = -float(max_py[i]) 
                 pos_weight[dpindex[i]] = pos_weight[dpindex[i]]*exp(v)
          
         ###### negative samples #####
           for i in range(num_of_samples):
                 v = float(max_py[i+num_of_samples]) 
                 neg_weight[dpindex[i]] = neg_weight[dpindex[i]]*exp(v)
          
         
           pos_w_sum = 0.0
           neg_w_sum = 0.0

           for i in range(pos_samples):
                 pos_w_sum += pos_weight[i]
  
           for i in range(neg_samples):
                 neg_w_sum += neg_weight[i]
           
           f = open('weight','w')

           for i in range(pos_samples):
                 pos_weight[i] = pos_weight[i]/pos_w_sum
                 tpw = str(pos_weight[i]) + '\n'
                 f.write(tpw) 

           for i in range(neg_samples):
                 neg_weight[i] = neg_weight[i]/neg_w_sum
                 tnw = str(neg_weight[i]) + '\n'
                 f.write(tnw)

           f.close()   
         
         ##################################################         
         ###### Perform AUC check on the Entire Set #######
         ##################################################
   
         param = ""
         m_train_y , m_train_x = myread(sep_data[max_index])

         print m_train_y[0]
         print m_train_x[0]
 
         m_test_y , m_test_x = myread(sep_data[max_index])  
         m_deci,m_model,m_py = get_pos_deci(m_train_y, m_train_x, m_test_y, m_test_x, param)
         print m_py
         eauc = plot_roc(m_deci,m_test_y)       
         print "AUC on entire set with window index = ",max_index," is : ",eauc
         svm_save_model('my.model',m_model)


         if max_index not in optimum:
                optimum.append(max_index)
         
         f = open('optimum','w')
         for m_i in optimum:
                t_s_i = str(m_i) + '\n'  
		f.write(t_s_i)

         f.close()
 
main()








